function ses_nstats()
% basic template file for data analysis routines.
% gets the current datac window selection, and does the analysis
% for all selected record groups/blocks
% automatically refreshes the data arrays from the disk file if necessary
% and also shows how to generate a figure.
% 8/23/05 P. Manis

% the main routine is just a wrapper to handle the multiple selection and
% catch execution errors.
%
% Note: when testing, set the testflag to 1 so that we DO NOT do try/catch.
% you need to do this to catch the errors in the _go routine.
%
global CONTROL
sf = getmainselection; % read the list of selected traces
testflag = 0;
if(sf > 0)
    pflag = getplotflag; % read teh plot flag (used in the routine template2)

    QueMessage('SES Noise Stats analysis routine', 1); % clear the que
    fprintf(1, 'File  Records   I (mean) I (std)  CV(I)   rate slope(sp/min)  delta-spikes\n');
    for i = 1:length(sf)
        if(testflag == 0)
            try
                template_go(sf(i), pflag);
            catch
                fprintf(1, '******ses_nstats Analysis Failed execution for sf = %d\n  ***[File: %s, records: %d - %d ]\n', ...
                    sf(i), CONTROL(sf(i)).filename, CONTROL(sf(i)).recbeg, CONTROL(sf(i)).recend);
            end;
        else
            template_go(sf(i), 0); % also turn off plot flag when testing
        end;

    end;
end;



%--------------The real analysis routine------------------

function template_go(sf, pflag)

% access to all of the data and parameters
global DFILE CONTROL VOLTAGE CURRENT

% set some control parameters
isrc = 1;
t0 = 0;
tmax = 1000;
spike_thresh = 0;

% this next routine reads all of the data and information we need,
% including calculating some standard stuff that is put into DP.
% you will want DP.time especially - it is the time base that goes
% with the data.
% 
% [DFILE, DP, err] = ana_setup(DFILE, sf);
% if(err ~= 0)
%     return;
% end;

%if(~isempty(ALLCH))
%    VOLTAGE = ALLCH{isrc};
%    CURRENT = ALLCH{isrc+1};
%end;
[DFILE] = synch_file(DFILE, sf);
if(isempty(DFILE) | isempty(VOLTAGE))
   fprintf('? could not get dfile');
   return;
end;

protocol=lower(CONTROL(sf).protocol);
QueMessage(sprintf('Processing %s, %s for ses_nstats', DFILE.filename, protocol));
RATES = DFILE.rate.*DFILE.nr_channel/1000;

[records,pts]=size(CURRENT);
time=make_time(DFILE);
tmax=max(max(time'));
if(DFILE.mode >= 5)
   wz=DFILE.ztime;
   w=find(diff(wz) < -12*60*60); % correct for possibility that someone actually records the data across midnight rollover.... (yes, it happened. 5/16/01 with Huijie's data.)
   if(~isempty(w))
      wz(w+1:end)=wz(w+1:end)+24*60*60;
   end;
   ZT = (wz-wz(1))/(60);
   cond_baseline = 5; % 5 min baseline
elseif length(DFILE.ztime) == length(RL)
   ZT = (DFILE.ztime-DFILE.ztime(1))/(60*1000);
   cond_baseline = 3;
else
   ZT = ones(length(RL), 1)*DFILE.cycle;
   cond_baseline = 5;
end;

% We need to find the spikes.
t0 = floor(10/RATES(1));
t1 = floor(900/RATES(1)); % time windows...

[first_spike, first_isi, nr_spikes, splat]=find_spikes2(DFILE, VOLTAGE, 100, 900, CONTROL(sf).thresh);

%now get the mean current and it's standard deviation.
t0 = floor(100/RATES(1));
t1 = floor(900/RATES(1)); % time windows...

tbase = floor(1/RATES(1));
for i = 1:size(VOLTAGE,1)
    ih(i) = mean(VOLTAGE(i,1:tbase)); % baseline voltage....
    imean = mean(CURRENT(i,t0:t1));
    istd = std(CURRENT(i,t0:t1));
    if(imean > 0)
        cv(i) = istd/imean;
    else
        cv(i) = NaN; % artifically set to 0 if it is not defined.
    end;
end;
DP.ZT = ZT;

[m, b, r, p] = linreg(DP.ZT', nr_spikes); % regression on the spike count for stability
deltas = (mean(nr_spikes(end-5:end) - mean(nr_spikes(1:5))))/mean(nr_spikes);

fprintf(1, '%s  %s   %7.2f  %7.2f  %8.3f  %8.2f  %8.3f\n', ...
    DFILE.filename, CONTROL(sf).reclist, imean, istd, mean(cv), m, deltas);

% set up a figure window (or you could grab an existing one)
hf = findobj('tag', 'Template'); % find the tagged window with it's name
if(isempty(hf)) % doesn't exist
    hf = figure();
    set(hf, 'tag', 'Template');
    set(hf, 'NumberTitle', 'off');
end;
figure(hf);


clf;% clear it.


% Define a text area and print some standard junk.
subplot('position', [0.1, 0.90, 0.5, 0.095]);
axis([0,1,0,1])
axis('off')
ht(1)=text(0,0.80,sprintf('%-12s R[%d:%d]     %-8s  [%s]',DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol, date), 'Fontsize', 10);
set(ht(1), 'interpreter', 'none'); % un-TeX the line - this is a filename and won't have tex chars, but might have an underscore.
text(0,0.6,sprintf('Solution:%-12s  gain:%4.1f  LPF:%4.1f kHz', CONTROL(sf).solution, DFILE.igain, DFILE.low_pass(1)), 'FontSize', 8);
%text(0,0.4,sprintf('Ihold:%6.2f %s    RMP: %6.2f %s, Rin: %8.3f M\\Omega', ...
%    CONTROL(sf).iHold,CONTROL(sf).I_Unit, CONTROL(sf).Rmp, CONTROL(sf).V_Unit, CONTROL(sf).Rin), 'FontSize', 8);


% define some subplots and plot some data

subplot('position', [0.1, 0.07, 0.45, 0.17]);
plot(DP.ZT, ih);
h1 = gca; % get a handle to the axes.
% label the plot and stuff if you need
xlabel('Time (min)');
ylabel('RMP (mV)');
box off % turn off the outer right and top axes marks.

subplot('position', [0.1, 0.27, 0.45, 0.17]);
plot(DP.ZT, nr_spikes);
hold on;
h1 = gca;
%xlabel('Time (ms)');
ylabel('Spike count');
y = m*DP.ZT+b;
plot(DP.ZT, y, 'r-');
box off

subplot('position', [0.1, 0.47, 0.45, 0.17]);
plot(DP.ZT, first_spike);
h1 = gca;
%xlabel('Time (min)');
ylabel('FSL');
box off

subplot('position', [0.1, 0.67, 0.45, 0.17]);
plot(DP.ZT, cv);
set(gca, 'YLim', [0 1]);
h1 = gca;
%xlabel('Time (min)');
ylabel('CV');
box off

subplot('position', [0.66, 0.47, 0.32, 0.40]);
ispec(gca, [100 900], 1);
xlabel('F (Hz)');
ylabel('P (sqrt-Hz)');

if(pflag) % plotflag and enable figure printing
    orient landscape % make the graph landscape on the paper.
    print;
end;



